Filtered model output statistics (FMOS)

ABSTRACT

A computer-implemented method to provide stabilized and spatially smooth regression coefficients for weather forecast error correction from small training data sets. In accordance with the present invention, an MOS estimate of the regression coefficient calculated from a small data set can be optimally combined with a smooth prior estimate of the regression coefficient, an estimate of the spatial error covariance of that prior estimate, and an estimate of the spatial error covariance of the MOS estimate. The result is a filtered MOS (FMOS) regression coefficient which can be used to more accurately estimate and correct errors in weather forecasts even using only small data sets.

CROSS-REFERENCE

This application claims the benefit of priority based on U.S.Provisional Patent Application No. 61/181,686 filed on May 28, 2009, theentirety of which is hereby incorporated by reference into the presentapplication.

TECHNICAL FIELD

The present invention relates to weather forecasting and statisticalcorrections of forecast errors.

BACKGROUND

Weather forecasts are of high value to a wide range of activities inagriculture, energy, transportation, defense and supply chainmanagement. See K. Rollins and J. Shaykewich, “Using willingness-to-payto assess the economic value of weather forecasts for multiplecommercial sectors,” Meteorol. Appl. Vol. 10, pp. 31-38 (2003); Y. Zhu,Z. Toth, R. Wobus, D. Richardson, and K, Mylne, “The Economic Value ofEnsemble-Based Weather Forecasts,” BAMS, January 2003, pp. 73-83. Ofparticular value are weather forecasts tailored for the particularweather related decisions that individuals, companies and organizationsmust make.

Such tailored weather forecasts are often provided by private weathercompanies such as those listed on the government web site“www(dot)weather(dot)gov/im/more.htm”.

The starting point for these weather forecasts is a computer simulationor Numerical Weather Prediction (NWP) in which approximations to theequations governing the evolution of the atmosphere propagate a recentlyestimated weather state forward in time. In general, substantialimprovements to forecasts made in this way can be achieved by comparinghistories of forecasts with corresponding histories of verifyingobservations. By mathematically modeling the differences between the NWPand observations, it has been shown that one can produce a forecast thatis significantly better than that from NWP.

Glahn and Lowry describe how the theory of multivariate linearregression can be used to correct systematically predictable aspects ofweather forecast models such as the NWP, and refer to their weatherforecast error corrector as Model Output Statistics (MOS). See H. R.Glahn and D. A. Lowry, “The Use of Model Output Statistics (MOS) inObjective Weather Forecasting,” J. Appl. Meteor., Vol. 11, Issue 8, pp.1203-1211 (1972). The MOS forecast error corrector is estimated from ahistorical record of forecasts together with corresponding verifyingobservations. Wilson and Vallée (2002) suggest that the historicalrecord needs to contain at least two years worth of daily weatherforecasts. L. J. Wilson and M. Vallée, “The Canadian Updateable ModelOutput Statistics (UMOS) System: Design and Development Tests,” Wea.Forecasting, Vol. 17, Issue 2, pp. 206-222 (2002).

MOS methods assume that on the kth day or event, the ith model forecastvariable f_(i) ^(k) is related to the corresponding verifying analysisor observation variable y_(i) ^(k) via the stochastic relationf _(i) ^(k) =a _(i) y _(i) ^(k) +b _(i)+ε_(i) ^(k)  (1)where the true regression coefficients a_(i) and b_(i) of the modelforecast variable f_(i) ^(k) can be considered to be the “slope”parameter and the “intercept” parameter, respectively, of a plot of therelation between f and y, and ε_(i) ^(k) is a random number associatedwith the kth forecast that is statistically independent of y_(i) ^(k).

With a finite sample of k=1, M realizations, using MOS methods one canobtain estimates a_(i) ^(MOS) and b_(i) ^(MOS) of the regressioncoefficients a_(i) and b_(i) using the relations

$\begin{matrix}{{a_{i}^{MOS} = \frac{\left\lbrack {\sum\limits_{k = 1}^{M}\;{\left( {f_{i}^{k} - {\overset{\_}{f}}_{i}} \right)\left( {y_{i}^{k} - {\overset{\_}{y}}_{i}} \right)}} \right\rbrack}{\sum\limits_{k = 1}^{M}\;{\left( {y_{i}^{k} - {\overset{\_}{y}}_{i}} \right)\left( {y_{i}^{k} - {\overset{\_}{y}}_{i}} \right)}}},{b_{i}^{MOS} = {{\overset{\_}{f}}_{i} - {a_{i}^{MOS}{\overset{\_}{y}}_{i}}}},{{{where}\mspace{14mu}{\overset{\_}{f}}_{i}} = {\frac{1}{M}{\sum\limits_{k = 1}^{M}f_{i}^{k}}}},{{\overset{\_}{y}}_{i} = {\frac{1}{M}{\sum\limits_{k = 1}^{M}y_{i}^{k}}}}} & (2)\end{matrix}$where f_(i) ^(k) is the ith model forecast variable and y_(i) ^(k) isthe verifying analysis/observation variable in Equation (1) above and f_(i) and y _(i) are the average values over M realizations.

If two years of historical data are available for the forecast model ofinterest, MOS can be relied on to provide significant forecastimprovements. However, in the continuing effort to improve weatherforecasting models, significant changes are frequently made to theforecasting model. The MOS equations developed for an old model may beentirely inappropriate for the changed new model. Consequently, afterevery model change, new MOS equations need to be developed. For example,typical model changes have included increases in model resolution,changes in the representation of sub sub-grid scale physics, and changesin the forcings associated with radiation.

If sufficient manpower and computer resources were available, ahistorical record of the performance of the new system required for newMOS equations could be generated fairly quickly by running the newsystem on historical data. However, such resources are generallyunavailable to the responsible parties. Instead, a historical record oftwo years is usually obtained by archiving forecasts and thecorresponding verifying analyses and forecasts in real time. Using sucha method, a two-year historical record would require two years' time tocompile.

This state of affairs led Wilson and Vallée (2002), supra, to propose anUpdateable Model Output Statistics (UMOS) system that blended regressionequations based on earlier versions of the model with equations based onmore recent versions of the model. UMOS implicitly assumes thatsignificant model changes are typically separated by more than twoyears. However, this assumption is not satisfied at many forecastingcenters. For example, the limited area weather forecasting models of theU.S. Navy are typically deployed in a specific theatre for a periodbetween two weeks and two years. In principle, new MOS equations wouldneed to be developed for each unique theatre in which the forecastingmodel is deployed. Hence, the Navy's deployment periods are too shortfor the estimated regression coefficients to stabilize, and any attemptto estimate the regression coefficients using MOS or UMOS with such asmall data set will result in very noisy regression coefficients.

One alternative would be to assume that regression coefficients wereidentical in various sub-regions of the model. This was assumed by, forexample, M. J. Schmeits, K. J. Kok, and D. H. P. Vogelezang,“Probabilistic Forecasting of (Severe) Thunderstorms in the NetherlandsUsing Model Output Statistics.” Wea. Forecasting, Vol. 20, Issue 2, pp.134-148 (2005). Although this approach will give stable coefficients forsmall data training sets it does not converge to the (superior) MOScoefficients in the limit of infinite training data. In addition, theMOS equations are likely to have erroneous jumps between sub-regions.

Thus, there is a need for a method for producing less noisy MOSregression coefficients that can be used with small data sets.

SUMMARY

This summary is intended to introduce, in simplified form, a selectionof concepts that are further described in the Detailed Description. Thissummary is not intended to identify key or essential features of theclaimed subject matter, nor is it intended to be used as an aid indetermining the scope of the claimed subject matter. Instead, it ismerely presented as a brief overview of the subject matter described andclaimed herein.

The present invention provides a computer-implemented method to providestabilized and spatially smooth regression coefficients for weatherforecast error correction from small training data sets. In accordancewith the present invention, the MOS estimates a_(i) ^(MOS) and b_(i)^(MOS) of regression coefficients calculated from a small data set canbe optimally combined with a smooth prior estimate of the regressioncoefficient, an estimate of the spatial error covariance of that priorestimate, and an estimate of the spatial error covariance of the MOSestimate to produce a set of Filtered MOS (FMOS) regression coefficientsa^(FMOS) and b^(FMOS):a ^(FMOS) =a ^(prior) =G _(a)(G _(a) =R _(a))⁻¹(a ^(MOS) −a ^(prior))b ^(FMOS) =b ^(prior) =G _(b)(G _(b) =R _(b))⁻¹(b ^(MOS) −b ^(prior))where(a ^(MOS))^(T) =[a ₁ ^(MOS) ,a ₂ ^(MOS) , . . . ,a _(n) ^(MOS)]and(b ^(MOS))^(T) =[b ₁ ^(MOS) ,b ₂ ^(MOS) , . . . ,b _(n) ^(MOS)]are n-vectors listing all the estimates of a₁ ^(MOS) and b_(i) ^(MOS)for i=1, . . . , n, n being the number of grid points defining the fieldof interest a^(prior) and b^(prior) are n-vectors listing priorestimates or guesses of the regression coefficients a and b over thatfield, G_(a) and G_(b) define the error covariances of a^(prior) andb^(prior) and R_(a) and R_(b) define the error covariance of a^(MOS) andb^(MOS).

These filtered regression coefficients a^(FMOS) and b^(FMOS) obtained inaccordance with the present invention can be used to more accuratelyestimate and correct errors in weather forecasts even using only smalldata sets.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 depicts an exemplary random random meridional wind field (v)generated using simulated observed idealized weather fields.

FIGS. 2A-2H illustrate aspects of the differences between the MOS andFMOS estimates of the slope coefficient a obtained from 16trials/events. FIGS. 2A and 2B illustrate an exemplary set of MOSestimates (prior art) of the true a-field and the error in thisestimate, respectively. FIG. 2C shows a plot of an exemplary truea-field. FIG. 2D shows the true spatial covariance function of thea-field. FIGS. 2E and 2F illustrate estimated error covariances of theprior and MOS estimates of a, respectively. FIGS. 2G and 2H,respectively, illustrate an exemplary set of Filtered MOS (FMOS)estimates of a obtained in accordance with the present invention and theerrors in the FMOS estimates as compared to the true a-field.

FIGS. 3A-3H illustrate corresponding aspects of the differences betweenthe MOS and the FMOS estimates of the slope coefficient b. Thus, FIGS.3A and 3B, illustrate an exemplary set of MOS estimates (prior art) ofthe true b-field and the error in this estimate, respectively. FIG. 3Cshows a plot of an exemplary true b-field. FIG. 3D shows the truespatial covariance function of the b-field. FIGS. 3E and 3F illustrateestimated error covariances of the prior and MOS estimates of b,respectively. FIGS. 3G and 3H, respectively, illustrate an exemplary setof FMOS estimates of b obtained in accordance with the present inventionand the errors in the FMOS estimates as compared to the true b-field.

DETAILED DESCRIPTION

The aspects and features of the present invention summarized above canbe embodied in various forms. The following description shows, by way ofillustration, combinations and configurations in which the aspects andfeatures can be put into practice. It is understood that the describedaspects, features, and/or embodiments are merely examples, and that oneskilled in the art may utilize other aspects, features, and/orembodiments or make structural and functional modifications withoutdeparting from the scope of the present disclosure.

The present invention provides a computer-implemented method forstatistically deriving weather forecast improvement functions from smalltraining data sets. As will be appreciated by one skilled in the art, amethod for providing stabilized and spatially smooth regressioncoefficients for weather forecast error correction in accordance withthe present invention can be accomplished by executing one or moresequences of instructions contained in computer-readable program coderead into a memory of one or more general or special-purpose computersconfigured to execute the instructions, wherein MOS regressioncoefficients that have large errors because of the paucity of trainingdata are transformed into a set of filtered modeled output statistics(FMOS) regression coefficients using information about the spatialcovariance of the prior estimate of the coefficients and the spatialerror covariance of the MOS estimate of the coefficients. The additionof the spatial covariance information causes the FMOS coefficients to bemuch more accurate than the noisy MOS coefficients from which they werederived.

As described in more detail below, the FMOS technique of the presentinvention stabilizes and spatially smoothes noisy regressioncoefficients produced by conventional MOS techniques using small datasets by optimally combining a data-poor MOS estimate with a priorestimate based on an assumption, for example, that model forecasts areunbiased. To optimally combine these two estimates, information aboutthe spatial covariance of the prior estimate of the coefficients and thespatial error covariance of the MOS estimate of the coefficients isused. An estimate of the spatial error covariance G of the priorestimate can be obtained from a spectral analysis of the (noisy)difference between the data-poor MOS estimate and the prior estimate,and an estimate of the spatial error covariance R of the data-poor MOSestimate can obtained from a spatial analysis of the equations governingthe grid-point wise error variance of the MOS equations. In accordancewith the present invention, these error covariances G and R are combinedwith the MOS coefficients to produce a set of filtered MOS regressioncoefficients which can be used to more accurately estimate and correcterrors in weather forecasts even using only small data sets.

As noted above, at the kth event (e.g., an observation day) a forecastf_(i) ^(k) at the ith grid point from the true state y_(i) ^(k) can beestimated from Equation (1) set forth above and replicated here forconvenience:f _(i) ^(k) =a _(i) y _(i) ^(k) +b _(i)+ε_(i) ^(k).  (1)where the regression coefficients a_(i) and b_(i) of the modeledforecast variable f_(i) ^(k) can be considered to be the “slope”parameter and the “intercept” parameter, respectively, of a plot of therelation between f and y, and ε_(i) ^(k) is a random number associatedwith the kth forecast that is statistically independent of y_(i) ^(k).Although Equation (1) expresses the forecast as function of theobservation rather than the observation as a function of the forecast,one skilled in the art would readily understand that the methodsdescribed herein can be applied equally well to counterparts of Equation(1) that express the observation as a function of the forecast. It maybe noted, however, that the form of Equation (1) presented herein may bepreferred when one is attempting to produce reliable probabilisticforecasts from ensemble forecasts. See C. H. Bishop and K. T. Shanley,“Bayesian Model Averaging's problematic treatment of extreme weather anda paradigm shift that fixes it,” Mon. Wea. Rev. 136, pp. 4641-4652(2008).

Also as noted above, Equation (1) contains a set of estimated regressioncoefficients a_(i) ^(MOS) and b_(i) ^(MOS) can be calculated usingEquation (2) previously presented above and replicated here forconvenience:

$\begin{matrix}{{a_{i}^{MOS} = \frac{\left\lbrack {\sum\limits_{k = 1}^{M}\;{\left( {f_{i}^{k} - {\overset{\_}{f}}_{i}} \right)\left( {y_{i}^{k} - {\overset{\_}{y}}_{i}} \right)}} \right\rbrack}{\sum\limits_{k = 1}^{M}\;{\left( {y_{i}^{k} - {\overset{\_}{y}}_{i}} \right)\left( {y_{i}^{k} - {\overset{\_}{y}}_{i}} \right)}}},{b_{i}^{MOS} = {{\overset{\_}{f}}_{i} - {a_{i}^{MOS}{\overset{\_}{y}}_{i}}}},{{{where}\mspace{14mu}{\overset{\_}{f}}_{i}} = {\frac{1}{M}{\sum\limits_{k = 1}^{M}f_{i}^{k}}}},{{\overset{\_}{y}}_{i} = {\frac{1}{M}{\sum\limits_{k = 1}^{M}y_{i}^{k}}}}} & (2)\end{matrix}$

However, there often are significant errors in the estimated regressioncoefficients a_(i) ^(MOS) and b_(i) ^(MOS), particularly in cases whereonly a small data set is present.

The present invention provides a method to improve on the MOS estimatesof the regression coefficients a_(i) and b_(i), which in turn can beused to provide more accurate weather forecasts and modeling, even incases where only small data sets are available. In accordance with thepresent invention, each MOS estimated coefficient can be combined with asmooth prior estimate of that regression coefficient, an estimate of thespatial error covariance of that prior estimate, and an estimate of thespatial error covariance of the MOS estimate.

The result is a set of filtered MOS (FMOS) regression coefficientsa^(FMOS) and b^(FMOS):a ^(FMOS) =a ^(prior) =G _(a)(G _(a) =R _(a))⁻¹(a ^(MOS) −a ^(prior))b ^(FMOS) =b ^(prior) =G _(b)(G _(b) =R _(b))⁻¹(b ^(MOS) −b^(prior))  (3)where(a ^(MOS))^(T) =[a ₁ ^(MOS) ,a ₂ ^(MOS), . . . ,a _(n) ^(MOS)]and(b ^(MOS))^(T) =[b ₁ ^(MOS) ,b ₂ ^(MOS), . . . ,b _(n) ^(MOS)]are n-vectors listing all the estimates of a_(i) ^(MOS) and b_(i) ^(MOS)for i=1, . . . , n, n being the number of grid points defining the fieldof interest a^(prior) and b^(prior) are n-vectors listing priorestimates or guesses of the regression coefficients a and b over overthat field, G_(a) and G_(b) define the error covariances of a^(prior)and b^(prior) and R_(a) and R_(b) define the error covariance of a^(MOS)and b^(prior) as described in more detail below. The improved regressioncoefficients a^(FMOS) and b^(FMOS) can be used to more accuratelyestimate and correct errors in weather forecasts even using only smalldata sets.

A method for estimating a^(FMOS) and b^(FMOS) in accordance with thepresent invention will now be described in detail.

As described above, the error covariances G_(a) and G_(b) give the errorcovariance of the prior estimates a^(prior) and b^(prior). n accordancewith the invention, G_(a) and G_(b) can be defined asG _(a)=

(a ^(prior) −a)(a ^(prior) −a)^(T)

,G _(b)=

(b ^(prior) −b)(b ^(prior) −b)^(T)

,  (4)where a and b are n-vectors listing the true values of the coefficientsa_(i) and b_(i) and a^(prior) and b^(prior) are n-vectors listing priorestimates or guesses of the regression coefficients. As noted above,after a change in the weather forecast model, one would typically assumethat (a^(prior))^(T)=[1, 1, . . . , 1] and (b^(prior))^(T)=[0, 0, . . ., 0] because, in large measure, the purpose of model changes is to makeforecasts unbiased and to make the approximation f_(i) ^(k)≈y_(i) ^(k)as accurate as possible. Thus, in the examples presented here, a^(prior)and b^(prior) will be based on the assumption that the model isunbiased, where all of the elements of a^(prior) are set equal to unityand all of the elements of b^(prior) are set equal to zero. However, ingeneral a^(prior) and b^(prior) can represent any prior guess of theregression coefficients a and b.

An exemplary method for estimating the covariances G_(a) and G_(b) inaccordance with the present invention will now be described. Since theestimation method is the same for both G_(a) and G_(b), only the methodfor estimating G_(a) will be described here, as one skilled in the artwould understand how to estimate G_(b) from the description herein. Inaddition, one skilled in the art would readily appreciate that othercovariance matrices G_(a) and G_(b) and methods for estimating thereofalso can be used, and any appropriate covariance matrices G_(a) andG_(b) and/or methods for estimating thereof can be used in estimatingregression coefficients a^(FMOS) and b^(FMOS) in accordance with thepresent invention.

Let an n×n matrix E list the real set of discrete orthonormalsinusoidal/cosinusoidal basis functions for the domain of interest andE^(T) be the transpose of E. These basis functions can be readilyaccessed via Fast Fourier Transform subroutine libraries known in theart. If the domain under consideration lies on a sphere, then a basis ofspherical harmonics should be used rather than sinusoids. If theassumption of periodicity is invalid on a finite domain then a puresinusoid or pure cosinusoid basis should be used.

It can be assumed thatG _(a) =E(diag(g))E ^(T)  (5)

The term diag(g) in Equation (5) is a diagonal matrix whose non-zeroelements are given byg _(k,l) =Bexp(−b ²(k ² +l ²)),g ^(T) =[g _(k,l) ,k=1, . . . nk,l=1, . . . nl]  (6)where the subscripts k and l refer to a sinusoidal basis function withwavenumber k in the x-direction and wavenumber l in the y-direction.

The broader the covariance functions associated with the covariancematrices G_(a) and G_(b), the smaller the historical data set requiredfor FMOS to yield accurate regression coefficients. This is because whenthese covariance functions are broad, Equation (6) can remove the noisefrom the MOS estimates of a and b by spatially averaging the noisy MOSestimates. The breadth of the covariance functions depends on thepositive definite parameter b used in deriving the diagonal matrixelements g_(k,l), with a larger value of b leading to broader covariancefunctions and a smaller value of b leading to tighter covariancefunctions. The parameter B used in defining the diagonal matrix elementsg_(k,l) determines the amplitude of the covariance functions, withlarger (smaller) B values giving larger (smaller) amplitudes.

The prime objective of the estimation procedure described here is tochoose B and b so that the vector (a^(prior)−a) would appear to be arandom draw of an n-vector from a normal distribution with covarianceG_(a). In other words, g should be chosen so that the assumption(a ^(prior) −a)=E[diag(g)]^(1/2) ξ,ξ˜N(0,I)  (7)where ξ is some random normal vector with mean 0 and covariance I, is asplausible as possible.

To do this, we begin by taking the transformα=E ^(T)(a ^(prior) −a ^(MOS)),d=α⊙α  (8)where the symbol ⊙ indicates the elementwise vector product. Thetransform yields the vector α of spectral coefficients defininga^(prior) and a^(MOS). Since α is real and d=α⊙α, there are no negativeelements in d. For a two-dimensional domain, each element d_(k,l) givesthe square of the coefficient of the sinusoid with wavenumber k in thex-direction and wavenumber l in the y-direction.

To make the assumption in Equation (7) plausible, the parameters B and bin in Equation (6) can be chosen to minimize the function

$\begin{matrix}{{J = {\frac{1}{2}{\sum\limits_{k = 1}^{n_{k}}\;{\sum\limits_{l = 1}^{n_{l}}\left( {\left( g_{k,l} \right) - \left( d_{k,l} \right)} \right)^{2}}}}}{{by}\mspace{14mu}{setting}}{\frac{\partial J}{\partial B} = {{\sum\limits_{k = 1}^{n_{k}}\;{\sum\limits_{l = 1}^{n_{l}}{\frac{\partial g_{k,l}}{\partial B}\left( {\left( g_{k,l} \right) - \left( d_{k,l} \right)} \right)}}} = 0}}{and}{\frac{\partial J}{\partial\left( b^{2} \right)} = {{\sum\limits_{k = 1}^{n_{k}}\;{\sum\limits_{l = 1}^{n_{l}}{\frac{\partial g_{k,l}}{\partial\left( b^{2} \right)}\left( {\left( {fg}_{k,l} \right) - \left( d_{k,l} \right)} \right)}}} = 0.}}} & (9)\end{matrix}$

We also note that from Equation (6)

$\begin{matrix}{{\left( g_{k,l} \right) = {B\;{\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}}}{{\frac{\partial g_{k,l}}{\partial B} = {\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}},{\frac{\partial g_{k,l}}{\partial\left( b^{2} \right)} = {{- \left( {k^{2} + l^{2}} \right)}B\;{\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}}}}} & (10)\end{matrix}$so that using the relations in Equation (10) in Equation (9) gives us

$\begin{matrix}{{\frac{\partial J}{\partial B} = {{\sum\limits_{k = 1}^{n_{k}}\;{\sum\limits_{l = 1}^{n_{l}}{{\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}\left( {{B\;{\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}} - \left( d_{k,l} \right)} \right)}}} = 0}}\mspace{79mu}{and}} & (11) \\{\frac{\partial J}{\partial\left( b^{2} \right)} = {{\sum\limits_{k = 1}^{n_{k}}\;{\sum\limits_{l = 1}^{n_{l}}{{- \begin{pmatrix}{k^{2} +} \\l^{2}\end{pmatrix}}B\;{\exp\left\lbrack {- {b^{2}\begin{pmatrix}{k^{2} +} \\l^{2}\end{pmatrix}}} \right\rbrack}\begin{pmatrix}{{B\;{\exp\left\lbrack {{- b^{2}}\left( {k^{2} + l^{2}} \right)} \right\rbrack}} -} \\\left( d_{k,l} \right)\end{pmatrix}}}} = {\left. 0\Rightarrow{\sum\limits_{k = 1}^{n_{k}}\;{\sum\limits_{l = 1}^{n_{l}}{{- \begin{pmatrix}{k^{2} +} \\l^{2}\end{pmatrix}}{\exp\left\lbrack {- {b^{2}\begin{pmatrix}{k^{2} +} \\l^{2}\end{pmatrix}}} \right\rbrack}\begin{pmatrix}{{B\;{\exp\left\lbrack {{- b^{2}}\left( {k^{2} + l^{2}} \right)} \right\rbrack}} -} \\\left( d_{k,l} \right)\end{pmatrix}}}} \right. = 0}}} & (12)\end{matrix}$

The derivatives

$\frac{\partial J}{\partial B}\mspace{14mu}{and}\mspace{14mu}\frac{\partial J}{\partial\left( b^{2} \right)}$can then be used in any suitable minimization algorithm to find thevalues of B and b² that minimize J. With B and b defined in this way,the optimum error covariance G_(a)=E(diag(g))E^(T) can then be estimatedusing Equations (5) and (6) described above.

Similarly to G_(a), the error covariance matrix G_(b) is estimated asG_(b)=E(diag(g^(b)))E^(T), where E is an n×n matrix listing a real setof discrete orthonormal sinusoidal/cosinusoidal basis functions fordomain of interest, E^(T) is the transpose of E, diag(g^(b)) is adiagonal matrix whose non-zero elements are given byg ^(b) _(k,l) =Aexp(−a ²(k ² +l ²)),g ^(bT) =[g ^(b) _(k,l) ,k=1, . . . ,n _(k) ,l=1, . . . ,n _(l)]  (12.1)where k and l are sinusoidal basis functions with wavenumber k in thex-direction and wavenumber l in the y-direction, and where A and aminimize the function

$\begin{matrix}{{{J^{b} = {\frac{1}{2}{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}\left( {\left( g_{k,l}^{b} \right) - \left( d_{k,l}^{b} \right)} \right)^{2}}}}}{{by}\mspace{14mu}{setting}}\frac{\partial J^{b}}{\partial A} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{\frac{\partial g_{k,l}^{b}}{\partial A}\left( {\left( g_{k,l}^{b} \right) - \left( d_{k,l}^{b} \right)} \right)}}} = 0}}{and}{{\frac{\partial J^{b}}{\partial\left( a^{2} \right)} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{\frac{\partial g_{k,l}^{b}}{\partial\left( a^{2} \right)}\left( {\left( g_{k,l}^{b} \right) - \left( d_{k,l}^{b} \right)} \right)}}} = 0}},}} & (12.2)\end{matrix}$where d^(b) _(k,l) is the element of the vector d^(b) corresponding towavenumber k (x-direction) and wavenumber l (y-direction) whered ^(b) =[E ^(T)(b ^(prior) −b ^(MOS))]⊙[E ^(T)(b ^(prior) −b^(MOS))].  (12.2)

Note that the aforementioned derivatives of J^(b) take the specific form

$\begin{matrix}{{\frac{\partial J^{b}}{\partial A} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}\left( {{A^{b}{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}} - \left( d_{k,l}^{b} \right)} \right)}}} = 0}}\mspace{79mu}{and}} & (12.3) \\{\frac{\partial J^{b}}{\partial\left( a^{2} \right)} = {\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{{- \left( {k^{2} + l^{2}} \right)}A\;{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}\left( {{{A\;{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}} - \left. \quad\left( d_{k,l}^{b} \right) \right)} = {\left. 0\Rightarrow{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{{- \left( {k^{2} + l^{2}} \right)}{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}\left( {{A\;{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}} - \left( d_{k,l}^{b} \right)} \right)}}} \right. = 0.}} \right.}}}} & (12.4)\end{matrix}$

The derivatives

$\frac{\partial J^{b}}{\partial A}\mspace{14mu}{and}\mspace{14mu}\frac{\partial J^{b}}{\partial\left( a^{2} \right)}$can then be used in any suitable minimization algorithm to find thevalues of A and a² that minimize J^(b). With A and a defined in thisway, the vector g^(b) and the corresponding optimum error covarianceG_(b)=E(diag(g^(b)))E^(T) are also defined.

As can be seen from Equation (3) and as described above, in addition toG_(a) and G_(b), the method for estimating improved forecast errorregression coefficients a^(FMOS) and b^(FMOS) in accordance with thepresent invention also utilizes a second set of error covarianceequations R_(a) and R_(b), which define the error covariance of the MOSestimates of the regression coefficients a^(MOS) and b^(MOS):R _(a)=

(a ^(MOS) −a)(a ^(MOS) −a)^(T)

,R _(b)=

(b ^(MOS) −b)(b ^(MOS) −b)^(T)

  (13)

It should be noted that unlike the error covariances G_(a) and G_(b),each of the error covariance matrices R_(a) and R_(b) must be estimatedseparately.

As described in more detail below, these spatial error covariancematrices can be estimated from an analysis of the known equationsdefining the error variance of estimates of the regression coefficientsat each grid point. The aspects of this analysis relating to spatialcovariances are new and are not found in the prior art. In addition, themagnitude of the elements in R_(a) and R_(b) tend to zero as the lengthof the historical record M used as a training data set tends to infinityso that in this limit, b^(FMOS)=b^(FMOS)=b and a^(FMOS)=a^(MOS)=a.

An exemplary method for estimating R_(a) and R_(b) is presented below.However, as in the case of G_(a) and G_(b), other forms of R_(a) andR_(b) and methods for estimating thereof may exist and any appropriateform and/or method for estimation may be used in the method forestimating the regression coefficients a^(FMOS) and b^(FMOS) inaccordance with the present invention.

We begin with the estimation of R_(b). From Equations (1) and (3) setforth above,f _(i) =a _(i) y _(i) +b _(i)+ ζ _(i) =>b _(i) = f _(i) −a _(i) y _(i) −ζ _(i).  (14)since b_(i) ^(MOS)= f _(i)−a_(i) y _(i), b_(i) ^(MOS)−b_(i)= ζ _(i),where the overbar indicates a sample mean over the M training events,e.g.,

${\overset{\_}{f}}_{i} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}{f_{i}^{m}.}}}$

In conventional methods, the element {R_(b)}_(ij) on the ith row and jthcolumn of R_(b) can be estimated using the values of b_(i) ^(MOS) and(b_(i) ^(MOS)−b_(i)) from Equation (14) above:

$\begin{matrix}\begin{matrix}{\left\{ R_{b} \right\}_{ij} = \left\langle {\left( {b_{i}^{MOS} - b_{i}} \right)\left( {b_{j}^{MOS} - b_{j}} \right)} \right\rangle} \\{= \left\langle {{\overset{\_}{\zeta}}_{i}{\overset{\_}{\zeta}}_{j}} \right\rangle} \\{= \frac{\left\langle {\zeta_{i}\zeta_{j}} \right\rangle}{M}} \\{\approx \frac{\frac{1}{M - 1}{\sum\limits_{m = 1}^{M}\;{\left( {\zeta_{i}^{m} - {\overset{\_}{\zeta}}_{i}} \right)\left( {\zeta_{j}^{m} - {\overset{\_}{\zeta}}_{j}} \right)}}}{M}}\end{matrix} & (15)\end{matrix}$where the triangular brackets indicate the expectation operator known inthe art.

Using the approximation ζhd i^(m)≈f_(i) ^(m)−a_(i) ^(MOS)y_(i)^(m)−b_(i) ^(MOS) in Equation (15), an estimate of the covariancebetween the error of the bias estimate at the ith grid point and thebias estimate at the jth grid point can be obtained.

The approximation given by Equation (15) is a trivial extension offormulae well known in the field of multivariate regression. However,this approximation is extremely poor when the number of data sets M issmaller than the number of grid-points in the horizontal domain, andunfortunately, this data-poor state of affairs is what would typicallybe expected in meteorological applications.

Consequently, in accordance with the present invention, R_(b) is assumedto be diagonal in spectral space and can be estimated using the relation

$\begin{matrix}{R_{b} \approx {E\left\{ {{diag}\left\{ {E^{T}\left\{ \frac{\frac{1}{M - 1}{\sum\limits_{m = 1}^{M}\;{\left( {\zeta^{m} - \overset{\_}{\zeta}} \right)\left( {\zeta^{m} - \overset{\_}{\zeta}} \right)}}}{M} \right\} E} \right\}} \right\} E^{T}}} & (16)\end{matrix}$where ξ^(m)≈f^(m)−a^(MOS)⊙y^(m)−b^(MOS) gives the vector of differencesbetween the forecast and the MOS prediction of the forecast from the mthevent and

$\underset{\_}{\overset{\_}{\zeta}} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}{{\underset{\_}{\zeta}}^{m}.}}}$The symbol “⊙” denotes the elementwise vector product. The symbol{R_(b)}_(ij) hereafter denotes the element of the matrix R_(b) lying onits ith row and jth column

Note that in both of Equations (15) and (16), the magnitude of theelements of R_(b) approaches zero as M tends to infinity; i.e., theerror in the MOS coefficients a^(MOS) and b^(MOS) as compared to thetrue regression coefficients a and b goes to zero as the number oftraining events M goes to infinity. However, when the number of trainingdata sets M is significantly smaller than the number of grid points ofinterest, this approximation of the regression coefficient R_(b)provides results that are far superior than that obtained usingconventional methods.

As noted above, R_(a) and R_(b) must be estimated separately, and sohaving obtained an estimate of R_(b), we can now estimate R_(a).

As described above, the forecast variable f_(i) ^(m) is related to thecorresponding verifying analysis/observation variable y_(i) ^(m) byf_(i) ^(m)=a_(i)y_(i) ^(m)+b_(i)+ζ_(i) ^(m). Therefore,(f _(i) −{overscore (f)} _(i))(y _(i) −{overscore (y)} _(i))=a _(i) (y_(i) −{overscore (y)} _(i))² + ζ_(i)(y _(i) −{overscore (y)} _(i))  (17)where the barred variables f _(i) and y _(i) are mean variables asdescribed above.

Rearranging Equation (17), we geta _(i)=[ (f _(i) −{overscore (f)} _(i))(y _(i) −{overscore (y)} _(i))−ζ_(i)(y _(i) −{overscore (y)} _(i))]/ (y _(i) −{overscore (y)} _(i))²  (18)

From Equation (2) described above, the value of a_(i) ^(MOS) can bederived:a _(i) ^(MOS)=[ (f _(i) −{overscore (f)} _(i))(y _(i) −{overscore (y)}_(i))]/ (y _(i) −{overscore (y)} _(i))² ,  (19)

Using Equations (18) and (19) above, we can derive a_(i) ^(MOS)−a_(i):a _(i) ^(MOS) −a _(i)= ζ_(i)(y _(i) −{overscore (y)} _(i))/ (y _(i)−{overscore (y)} _(i))²   (20)and therefore can derive R_(a):

$\begin{matrix}\begin{matrix}{\left\{ R_{a} \right\}_{ij} = \left\langle {\left( {a_{i}^{MOS} - a_{i}} \right)\left( {a_{j}^{MOS} - a_{j}} \right)} \right\rangle} \\{= \left\langle {\overset{\_}{\zeta_{i}\left( {y_{i} - {\overset{\_}{y}}_{i}} \right)}{\overset{\_}{\zeta_{j}\left( {y_{j} - {\overset{\_}{y}}_{j}} \right)}/\left\lbrack {\overset{\_}{\left( {y_{i} - {\overset{\_}{y}}_{i}} \right)^{2}}\overset{\_}{\left( {y_{j} - {\overset{\_}{y}}_{j}} \right)^{2}}} \right\rbrack}} \right\rangle} \\{= \left\langle \frac{\left\{ {\left\lbrack {\frac{1}{M}{\sum\limits_{m = 1}^{M}{\zeta_{i}^{m}\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)}}} \right\rbrack\left\lbrack {\frac{1}{M}{\sum\limits_{k = 1}^{M}{\zeta_{j}^{k}\left( {y_{j}^{k} - {\overset{\_}{y}}_{j}} \right)}}} \right\rbrack} \right\}}{\left\lbrack {\frac{1}{M}{\sum\limits_{m = 1}^{M}\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)^{2}}} \right\rbrack\left\lbrack {\frac{1}{M}{\sum\limits_{k = 1}^{M}\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)^{2}}} \right\rbrack} \right\rangle} \\{\approx \frac{\left\{ \left\lbrack {\frac{1}{M^{2}}{\sum\limits_{k = 1}^{M}{\sum\limits_{m = 1}^{M}\left\langle {\zeta_{i}^{m}{\zeta_{j}^{k}\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)}\left( {y_{j}^{k} - {\overset{\_}{y}}_{j}} \right)} \right\rangle}}} \right\rbrack \right\}}{\left\lbrack {\frac{1}{M}{\sum\limits_{m = 1}^{M}\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)^{2}}} \right\rbrack\left\lbrack {\frac{1}{M}{\sum\limits_{k = 1}^{M}\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)^{2}}} \right\rbrack}}\end{matrix} & (21)\end{matrix}$

Using the facts that

y_(i) ^(m)ζ_(i) ^(m)

=

y_(i) ^(m)ζ_(i) ^(k)

=0 and

ζ_(i) ^(m)ζ_(j) ^(k)

=δ_(mk)

ζ_(i) ^(m)ζ_(j) ^(m)

, Equation (21) simplifies to

$\begin{matrix}\begin{matrix}{\left\{ R_{a} \right\}_{ij} = \frac{\left\{ \left\lbrack {\frac{1}{M^{2}}{\sum\limits_{m = 1}^{M}\left( {\left\langle {\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)} \right\rangle\left\langle {\zeta_{i}^{m}\zeta_{j}^{m}} \right\rangle} \right)}} \right\rbrack \right\}}{\left\lbrack {\frac{1}{M}{\sum\limits_{m = 1}^{M}\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)^{2}}} \right\rbrack\left\lbrack {\frac{1}{M}{\sum\limits_{k = 1}^{M}\left( {y_{j}^{k} - {\overset{\_}{y}}_{j}} \right)^{2}}} \right\rbrack}} \\{= \frac{\left\{ \left\lbrack {\frac{1}{M}{\sum\limits_{m = 1}^{M}\left( {\left\langle {\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)} \right\rangle\frac{\left\langle {\zeta_{i}^{m}\zeta_{j}^{m}} \right\rangle}{M}} \right)}} \right\rbrack \right\}}{\left\lbrack {\frac{1}{M}{\sum\limits_{m = 1}^{M}\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)^{2}}} \right\rbrack\left\lbrack {\frac{1}{M}{\sum\limits_{m = 1}^{M}\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)^{2}}} \right\rbrack}}\end{matrix} & (22)\end{matrix}$

Using Equation (15) in Equation (22) gives

$\begin{matrix}\begin{matrix}{\left\{ R_{a} \right\}_{ij} = \frac{\left\langle {\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)} \right\rangle\left\{ R_{b} \right\}_{ij}}{\left\lbrack {\frac{1}{M}{\sum\limits_{m = 1}^{M}\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)^{2}}} \right\rbrack\left\lbrack {\frac{1}{M}{\sum\limits_{m = 1}^{M}\left( {y_{j}^{M} - {\overset{\_}{y}}_{j}} \right)^{2}}} \right\rbrack}} \\{\approx \frac{\left\{ {\sum\limits_{m = 1}^{M}\frac{\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)}{M - 1}} \right\}\left\{ R_{b} \right\}_{ij}}{\left\lbrack {\frac{1}{M - 1}{\sum\limits_{m = 1}^{M}\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)^{2}}} \right\rbrack\left\lbrack {\frac{1}{M - 1}{\sum\limits_{m = 1}^{M}\left( {y_{j}^{M} - {\overset{\_}{y}}_{j}} \right)^{2}}} \right\rbrack}}\end{matrix} & (23)\end{matrix}$

To get Equation (23) into matrix form, we first define the matrix Y suchthat

$\begin{matrix}{\left\{ Y \right\}_{ij} = \left\{ {\sum\limits_{m = 1}^{M}\frac{\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)}{M - 1}} \right\}} & (24)\end{matrix}$

Let the diagonal matrix D list the inverse of the diagonal elements of Ysuch thatD={diag(Y)}⁻¹  (25)

Equation (23) may then be rewritten in the matrix formR _(a) ≈D(Y⊙R _(b))D  (26)where ⊙ denotes the elementwise matrix product. Note that Equation (26)states that each element of R_(a) is directly proportional to acorresponding element of R_(b). Consequently, just as the elements ofR_(b) tend to zero as the number of training data events M tends toinfinity, the elements of R_(a) will also tend to zero as M→∞.

Thus, in accordance with the present invention, with the MOS regressioncoefficients, the prior estimates of the regression coefficients, andthe covariances G_(a) G_(b), R_(a), and R_(b) known, an improved set ofregression coefficients a_(i) and b_(i), i.e., a_(i) ^(FMOS) and b_(i)^(FMOS), can be obtained in accordance with Equation (3) set forthabove:a ^(FMOS) =a ^(prior) +G _(a)(G _(a) +R _(a))⁻¹(a ^(MOS) −a ^(prior)),b ^(FMOS) =b ^(prior) +G _(b)(G _(b) +R _(b))⁻¹(b ^(MOS) −b ^(prior))

These improved regression coefficients in turn can be used to provide animproved statistical relationship between the forecast at the ith gridpoint f_(i) ^(k) and the observation y_(i) ^(k) at the ith grid pointfrom the FMOS coefficients a_(i) ^(FMOS) and b_(i) ^(FMOS) at the ithgrid point where:f _(i) ^(k) =a _(i) ^(FMOS) y _(i) ^(k) +b _(i) ^(FMOS)+ε_(i)^(k).  (27)

Equation (27) represents an improved estimate of the relationshipbetween observations and forecasts than that obtained using MOScoefficients because FMOS coefficients are more accurate than MOScoefficients. Using the methods of Bishop and Shanley (2008), supra,Equation (27) can be used to improve the accuracy and reliability ofprobabilistic forecasts of weather.

It will be noted here that in some embodiments of the present invention,all of a^(MOS), b^(MOS), a^(prior), b^(prior), G_(a) and G_(b), andR_(a) and R_(b) may already have been computed and stored in an internalor external memory or other storage medium, so that they need simply tobe loaded as inputs into a computer programmed with appropriate softwarewhich can then estimate the a^(FMOS) and b^(FMOS) coefficients using theinput values. In other embodiments, the values of a^(MOS) and b^(MOS)and a^(prior) and b^(prior) may have been previously computed and loadedas inputs, with G_(a) and G_(b), R_(a) and R_(b) being computed as partof the process of estimating the a^(FMOS) and b^(FMOS) coefficients. Instill other embodiments, none of a^(MOS), b^(MOS), G_(a) and G_(b),R_(a) and R_(b) may have been already computed so that all of a^(MOS),b^(MOS), G_(a) and G_(b), R_(a) and R_(b) are computed during theprocess of estimating the a^(FMOS) and b^(FMOS) coefficients.

To test the FMOS approach, random correlated fields of true valuesy^(T)=[y₁y₂, . . . , y_(n)] were created. An example of such a randomfield of true values is shown in FIG. 1, which illustrates a simulatedtrue state of a 2-dimensional wind component in m/s across idealizedweather fields having a scale of 1000 km. From these true values,corresponding forecast values f^(T)=[f₁, f₂, . . . , f_(n)] were createdat each grid point using the true values of spatially varying fields ofa and b in Equation (1). To give these fields meteorological relevancethey were made to propagate from West to East and evolve in time in away that is qualitatively similar to the way meteorological fields of,for example, temperature, evolve with time. In this way, time series ofpairs of observations and forecasts that obey Equation (1) could beeasily created for testing the MOS and FMOS approaches for estimatingthe coefficients a and b from the data fields. The fields together withthe MOS fields were also used to test the aforementioned methods forestimating G_(a) and G_(b), R_(a), and R_(b).

The beneficial effects of using the method of the present invention areclearly illustrated by the plots shown in FIGS. 2A-2H and FIGS. 3A-3H.

As noted above, FIGS. 2A-2H illustrate aspects of the differencesbetween the MOS and FMOS estimates of the slope coefficient a obtainedfrom 16 trials/events. FIGS. 3A-3H illustrate corresponding aspectsbetween the MOS and FMOS estimates of the slope coefficient b.

Each of the 16 trials produces an observation y and a forecast f at eachgrid point of the model domain. From these 16 data pairs, attempts weremade to estimate the unknown coefficients using the MOS method accordingto Equation (2) set forth above. The results of these estimates areshown in FIGS. 2A-2B and 3A-3B.

FIGS. 2C and 3C are plots of an exemplary true a-field and an exemplarytrue b-field obtained using a minor adjustment of the method describedin the appendix of Bishop and Hodyss (2007). See C. H. Bishop and D.Hodyss, “Flow adaptive moderation of spurious ensemble correlations andits use in ensemble based data assimilation.” Quart. J. Roy. Met. Soc.133, pp. 2029-2044 (2007). FIGS. 2D and 3D shows the true spatialcovariance function of the a- and b-fields shown in FIGS. 2C and 3C,respectively, and are equivalent to an error covariance function of theprior estimates of a and b that sets a and b equal to unity at all gridpoints.

FIGS. 2E and 2F show plots of G_(a), the estimated error covariancefunction of the prior estimates of a (FIG. 2E), and R_(a), the estimatederror covariance function of the MOS estimates of a (FIG. 2F), withFIGS. 3E and 3F showing the corresponding plots for G_(b) and R_(b).

FIGS. 2A/3A and 2B/3B illustrate aspects of the MOS estimates of thecoefficients a and b obtained using the methods of the prior art. Notethe noisy nature of the a and b fields shown in FIGS. 2A and 3A ascompared with the true fields shown in FIGS. 2C and 3C. In contrast, theFMOS estimates of the coefficients a and b obtained using the methodsdescribed herein are smooth, and resemble the plots of the true a and bfields shown in FIGS. 2C and 3C.

In addition, a comparison of the error of the FMOS estimates of a and bshown in FIGS. 2H/3H with the error of the MOS estimates shown in FIGS.2B/3B shows that the FMOS errors are about an order of magnitude smallerthan the MOS estimate. A comparison of these Figures shows that themaximum errors incurred by the MOS estimate are more than 3-10 timeslarger than those incurred by the FMOS estimate. For example, the errorof the a_(i) ^(MOS) field illustrated in FIG. 2B shows error extremawith magnitudes greater than 0.2, and in fact shows mean square errors(mse) of the MOS estimates that are 19 times larger than the mse of theFMOS estimate shown in FIG. 2H. The MOS estimate of the regressioncoefficients based on 16 events is of little practical use because thenoise (error) is larger than the signal. On the other hand, the errorobtained from the FMOS estimates shown in FIG. 2F have maximummagnitudes of 0.04, which is much smaller than the 0.2 amplitudevariations of a_(i). Thus, as can readily be appreciated, the FMOSestimates of the regression coefficients in accordance with the presentinvention represent a significant improvement over the MOS estimates foruse with small data sets.

In addition, to test the statistical significance of the superiority ofFMOS over MOS, 16 independent trials were performed in which

ε_(i) ²

=10

(a_(i)y_(i))²

=10

(b_(i))²

. In other words, the variance of the random component of forecast errorε_(i) in Equation (27) was made to be 10 times larger than the varianceof the signal. It was found that MOS mse exceeded FMOS mse for the slopeparameter “a” by a factor between 3.5 and 57.5. The mean of thissuperiority factor over the 16 trials was 15.3. In the same set ofindependent trials, MOS mse exceeded FMOS mse for the interceptparameter “b” by a factor between 2.4 and 9.6. The mean of thissuperiority factor was 6.1. Note that FMOS was superior to MOS in everyone of the 16 trials for both parameters. The probability of thishappening by pure chance if FMOS was no better than MOS is0.5¹⁶=1.5×10⁻⁵. Hence, the null hypothesis that FMOS is not superior toMOS can be rejected with a confidence level of 99.999%.

An additional 16-trial test was performed in which the variance of therandom component of forecast error in Equation (27) was reduced by anorder of magnitude so that

ε_(i) ²

=

(a_(i)y_(i))²

=

(b_(i))²

. In this case, MOS mse exceeded FMOS mse for the slope parameter “a” bya factor between 3.8 and 16.7. The mean of this superiority factor was7.8. In the same set of independent trials, MOS mse exceeded FMOS msefor the intercept parameter “b” by a factor between 1.5 and 5.8. Themean of this superiority factor was 3.7. Note that FMOS was superior toMOS in every one of the 16 trials. Since the probability of thishappening by pure chance if FMOS was no better than MOS is0.5¹⁶=1.5×10⁻⁵. Hence, the null hypothesis that FMOS is not superior toMOS can again be rejected with a confidence level of 99.999%.

Comparison of the first 16 trials with the second set of 16 trials showsthat the extent of superiority of FMOS over MOS increases as themagnitude of the systematic part of forecast error is reduced relativeto the random part of forecast error.

Thus, the use of the FMOS method of the present invention has numerousadvantages. For example, its estimates of regression coefficients aremore accurate, particularly when the training data set is small.Moreover, if properly configured, the results using the FMOS method ofthe present invention will be equivalent to those using MOS in the limitof an infinite data training set. This is a desirable feature because,in the limit of a very large data set, MOS is guaranteed to recover thetrue coefficients.

It will be appreciated by one skilled in the art that one or moreaspects of a method for providing stabilized and spatially smoothregression coefficients for weather forecast error correction from smalltraining data sets as described herein can be accomplished by one ormore processors executing one or more sequences of one or morecomputer-readable instructions read into a memory of one or morecomputers from volatile or non-volatile computer-readable media capableof storing and/or transferring computer programs or computer-readableinstructions for execution by one or more computers. Volatile media caninclude a memory such as a dynamic memory in a computer. Non-volatilecomputer readable media that can be used can include a compact disk,hard disk, floppy disk, tape, magneto-optical disk, PROM (EPROM, EEPROM,flash EPROM), SRAM, SDRAM, or any other magnetic medium; punch card;paper tape; or any other physical medium such as a chemical orbiological medium.

Although particular embodiments, aspects, and features have beendescribed and illustrated, it should be noted that the inventiondescribed herein is not limited to only those embodiments, aspects, andfeatures. It should be readily appreciated that modifications may bemade by persons skilled in the art, and the present applicationcontemplates any and all modifications within the spirit and scope ofthe underlying invention described and claimed herein. Such embodimentsare also contemplated to be within the scope and spirit of the presentdisclosure.

What is claimed is:
 1. A computer-implemented method for estimating a set of regression coefficients for use in estimating a model forecast variable from a small data set, the method including the steps of: (a) receiving, using a computer, data of an n-vector a^(MOS) and an n-vector b^(MOS) the elements of a^(MOS) and b^(MOS) comprising a set of modeled output statistics (MOS) estimates of the true values of forecast correction coefficients a and b; (b) receiving, using the computer, data of an n-vector a^(prior) and an n-vector b^(prior), the elements of a^(prior) and b^(prior) comprising a set of prior estimates of forecast correction coefficients a and b; (c) estimating a spatial covariance matrix G_(a) for a^(prior) and a spatial covariance matrix G_(b) for b^(prior); (i) wherein the spatial covariance matrix G_(a) is estimated as G_(a)=E(diag(g))E^(T), where E is an n×n matrix listing a real set of discrete orthonormal basis functions for domain of interest, E^(T) is the transpose of E, and diag(g) is a diagonal matrix whose non-zero elements are given by g _(k,l) =Bexp(−b ²(k ² +l ²)) and g ^(T) =[g _(k,j) ,k=1, . . . ,n _(k) ,l=1, . . . ,n _(l)]  in which k and l are sinusoidal basis functions with a wavenumber k in the x-direction and a wavenumber l in the y-direction, and B and b minimize the function $J = {\frac{1}{2}{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}\left( {\left( g_{k,l} \right) - \left( d_{k,l} \right)} \right)^{2}}}}$ by  setting $\frac{\partial J}{\partial B} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{\frac{\partial g_{k,l}}{\partial B}\left( {\left( g_{k,l} \right) - \left( d_{k,l} \right)} \right)}}} = 0}$ and ${\frac{\partial J}{\partial\left( b^{2} \right)} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{\frac{\partial g_{k,l}}{\partial\left( b^{2} \right)}\left( {\left( g_{k,l} \right) - \left( d_{k,l} \right)} \right)}}} = 0}},$  where d_(k,l) is the element of the vector d=[E^(T)(a^(prior)−a^(MOS))]⊙[E^(T)(a^(prior)−a^(MOS))] corresponding to the wavenumber k and the wavenumber l, where the symbol ⊙ indicates the elementwise product, and where the derivatives of J are of the form $\mspace{79mu}{\frac{\partial J}{\partial B} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{{\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}\left( {{B\;{\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}} - \left( d_{k,l} \right)} \right)}}} = 0}}$      and ${\frac{\partial J}{\partial\left( b^{2} \right)} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{{- \left( {k^{2} + l^{2}} \right)}b\;{\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}\left( {{B\;{\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}} - \left( d_{k,l} \right)} \right)}}} = {\left. 0\Rightarrow{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{{- \left( {k^{2} + l^{2}} \right)}{\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}\left( {{B\;{\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}} - \left( d_{k,l} \right)} \right)}}} \right. = 0}}};$ and (ii) wherein the spatial covariance matrix G_(b) is estimated as G_(b)=E(diag(g^(b)))E^(T), where diag(g^(b)) is a diagonal matrix whose non-zero elements are given by g ^(b) _(k,l) =Aexp(−a ²(k ² +l ²)) and g ^(bT) =[g ^(b) _(k,l) k=1, . . . ,n _(k) ,l=1, . . . ,n _(l)]  where k and l are the sinusoidal basis functions with the wavenumber k and the wavenumber l, and where A and a minimize the function $J^{b} = {\frac{1}{2}{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}\left( {\left( g_{k,l}^{b} \right) - \left( d_{k,l}^{b} \right)} \right)^{2}}}}$ by  setting $\frac{\partial J^{b}}{\partial A} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{\frac{\partial g_{k,l}^{b}}{\partial A}\left( {\left( g_{k,l}^{b} \right) - \left( d_{k,l}^{b} \right)} \right)}}} = 0}$ and ${\frac{\partial J^{b}}{\partial\left( a^{2} \right)} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{\frac{\partial g_{k,l}^{b}}{\partial\left( a^{2} \right)}\left( {\left( g_{k,l}^{b} \right) - \left( d_{k,l}^{b} \right)} \right)}}} = 0}},$  where d^(b) _(k,l) is an element of the vector d^(b)=[E^(T)(b^(prior)−b^(MOS))]⊙[E^(T)(b^(prior)−b^(MOS))] corresponding to the wavenumber k and the wavenumber l and the derivatives of J^(b) take the form $\mspace{79mu}{\frac{\partial J^{b}}{\partial A} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}\left( {{A^{b}{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}} - \left( d_{k,l}^{b} \right)} \right)}}} = 0}}$      and $\frac{\partial J^{b}}{\partial\left( a^{2} \right)} = {\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{{- \left( {k^{2} + l^{2}} \right)}A\;{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}\left( {{{A\;{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}} - \left. \quad\left( d_{k,l}^{b} \right) \right)} = \left. 0\Rightarrow{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{{- \left( {k^{2} + l^{2}} \right)}{\exp\left\lbrack {{- \mspace{79mu}\left. \quad{a^{2}\left( {k^{2} + l^{2}} \right)} \right\rbrack}\left( {{{A\;\left. \quad{{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack} - \left( d_{k,l}^{b} \right)} \right)} = 0};} \right.} \right.}}}} \right.} \right.}}}$ (d) estimating, using the computer, data of a spatial error covariance matrix R_(a) for a^(MOS) and a spatial error covariance matrix R_(b) for b^(MOS), each element {R_(a)}_(ij) of R_(a) and each element {R_(b)}_(ij) of R_(b) being estimated over a set of m=1 to M events, M being smaller than a number of points of interest in a spatial grid, (i) wherein R_(b) is estimated using the relation $R_{b} \approx {E\left\{ {{diag}\left\{ {E^{T}\left\{ \frac{\frac{1}{M - 1}{\sum\limits_{m = 1}^{M}{\left( {{\underset{\_}{\zeta}}^{m} - \underset{\_}{\overset{\_}{\zeta}}} \right)\left( {{\underset{\_}{\zeta}}^{m} - \underset{\_}{\overset{\_}{\zeta}}} \right)}}}{M} \right\} E} \right\}} \right\} E^{T}}$ ξ^(m)≈f^(m)−a^(MOS)⊙y^(m)−b^(MOS) gives a vector of differences between a forecast [(f_(m))^(T)=[f₁ ^(m), f₂ ^(m), . . . , f_(n) ^(m)]] and a MOS prediction f_(m) ^(MOS)=a^(MOS)⊙y^(m)+b^(MOS) of a forecast from an mth event and ${\underset{\_}{\overset{\_}{\zeta}} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}{\underset{\_}{\zeta}}^{m}}}},$  and wherein {R_(b)}_(ij) denotes an element of the matrix R_(b) lying on its ith row and jth column; and (ii) wherein R_(a) is estimated as R_(a)≈D(Y⊙R_(b))D, where ${\left\{ Y \right\}_{ij} = \left\{ {\sum\limits_{m = 1}^{M}\frac{\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)}{M - 1}} \right\}},{{{and}\mspace{14mu} D} = \left\{ {{diag}(Y)} \right\}^{- 1}},$  and wherein each element {R_(a)}_(ij) of error covariance matrix R_(a) is estimated as $\begin{matrix} {\left\{ R_{a} \right\}_{ij} = \frac{\left\langle {\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)} \right\rangle\left\{ R_{b} \right\}_{ij}}{\left\lbrack {\frac{1}{M}{\sum\limits_{m = 1}^{M}\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)^{2}}} \right\rbrack\left\lbrack {\frac{1}{M}{\sum\limits_{m = 1}^{M}\left( {y_{j}^{M} - {\overset{\_}{y}}_{j}} \right)^{2}}} \right\rbrack}} \\ {{= \frac{\left\{ {\sum\limits_{m = 1}^{M}\frac{\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)}{M - 1}} \right\}\left\{ R_{b} \right\}_{ij}}{\left\lbrack {\frac{1}{M - 1}{\sum\limits_{m = 1}^{M}\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)^{2}}} \right\rbrack\left\lbrack {\frac{1}{M - 1}{\sum\limits_{m = 1}^{M}\left( {y_{j}^{M} - {\overset{\_}{y}}_{j}} \right)^{2}}} \right\rbrack}},} \end{matrix}$  where y _(i) is a sample mean of a verifying variable y_(i) ^(m) over a plurality of m events; and (e) transforming the data of a^(MOS), b^(MOS), a^(prior), b^(prior), G_(a), G_(b), R_(a), and R_(b) into a set of estimated regression coefficients a^(FMOS) and b^(FMOS), wherein a ^(FMOS) =a ^(prior) +G _(a)(G _(a) +R _(a))⁻¹(a ^(MOS) −a ^(prior)); and b ^(FMOS) =b ^(prior) +G _(b)(G _(b) +R _(b))⁻¹(b ^(MOS) −b ^(prior)).
 2. A computer-implemented method for estimating a set of regression coefficients for use in estimating a model forecast variable from a small data set, the method including the steps of: (a) receiving, using a computer, data of an n-vector a^(MOS) and an n-vector b^(MOS), the elements of a^(MOS) and b^(MOS) comprising a set of modeled output statistics (MOS) estimates of the true values of forecast correction coefficients a and b; (b) receiving, using the computer, data of an n-vector a^(prior) and an n-vector b^(prior), the elements of a^(prior) and b^(prior) comprising a set of prior estimates of forecast correction coefficients a and b; (c) receiving, using the computer, data of a spatial covariance matrix G_(a) for a^(prior) and a spatial covariance matrix G_(b) for b^(prior); (d) estimating, using the computer, data of a spatial error covariance matrix R_(a) for a^(MOS) and a spatial error covariance matrix R_(b) for b^(MOS) each element {R_(a)}_(ij) of R_(a) and each element {R_(b)}_(ij) of R_(b) being estimated over a set of m=1 to M events, M being smaller than a number of points of interest in a spatial grid, (i) wherein R_(b) is estimated using the relation $R_{b} \approx {E\left\{ {{diag}\left\{ {E^{T}\left\{ \frac{\frac{1}{M - 1}{\sum\limits_{m = 1}^{M}{\left( {{\underset{\_}{\zeta}}^{m} - \underset{\_}{\overset{\_}{\zeta}}} \right)\left( {{\underset{\_}{\zeta}}^{m} - \underset{\_}{\overset{\_}{\zeta}}} \right)}}}{M} \right\} E} \right\}} \right\} E^{T}}$ where ξ^(m)≈f^(m)−a^(MOS)⊙y^(m)−b^(MOS) gives a vector of differences between a forecast [(f_(m))^(T)=[f₁ ^(m), f₂ ^(m), . . . , f_(n) ^(m)]], where the symbol ⊙ indicates the elementwise product, and a MOS prediction f_(m) ^(MOS)=a^(MOS)⊙y^(m)+b^(MOS) of a forecast from an mth event and ${\underset{\_}{\overset{\_}{\zeta}} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}{\underset{\_}{\zeta}}^{m}}}},$  wherein {R_(b)}_(ij) denotes an element of the matrix R_(b) lying on its ith row and jth column; and (ii) wherein R_(a) is estimated as R_(a)≈D(Y⊙R_(b))D, where ${\left\{ Y \right\}_{ij} = \left\{ {\sum\limits_{m = 1}^{M}\frac{\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)}{M - 1}} \right\}},{and}$ D = {diag(Y)}⁻¹,  and wherein each element {R_(a)}_(ij) of error covariance matrix R_(a) is estimated as $\begin{matrix} {\left\{ R_{a} \right\}_{ij} = \frac{\left\langle {\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)} \right\rangle\left\{ R_{b} \right\}_{ij}}{\left\lbrack {\frac{1}{M}{\sum\limits_{m = 1}^{M}\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)^{2}}} \right\rbrack\left\lbrack {\frac{1}{M}{\sum\limits_{m = 1}^{M}\left( {y_{j}^{M} - {\overset{\_}{y}}_{j}} \right)^{2}}} \right\rbrack}} \\ {{\approx \frac{\left\{ {\sum\limits_{m = 1}^{M}\frac{\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)}{M - 1}} \right\}\left\{ R_{b} \right\}_{ij}}{\left\lbrack {\frac{1}{M - 1}{\sum\limits_{m = 1}^{M}\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)^{2}}} \right\rbrack\left\lbrack {\frac{1}{M - 1}{\sum\limits_{m = 1}^{M}\left( {y_{j}^{M} - {\overset{\_}{y}}_{j}} \right)^{2}}} \right\rbrack}},} \end{matrix}$  where y _(i) is a sample mean of a verifying variable y_(i) ^(m) over a plurality of m events; and (e) transforming data of a^(MOS), b^(MOS), a^(prior), b^(prior), G_(a), G_(b), R_(a), and R_(b) into a set of estimated regression coefficients a^(FMOS) and b^(FMOS), wherein a ^(FMOS) =a ^(prior) +G _(a)(G _(a) +R _(a))⁻¹(a ^(MOS) −a ^(prior)); and b ^(FMOS) =b ^(prior) +G _(b)(G _(b) −R _(b))⁻¹(b ^(MOS) −b ^(prior)).
 3. A computer-implemented method for estimating a set of regression coefficients for use in estimating a model forecast variable f_(i) from a small data set, including the steps of: (a) receiving, using a computer, data of an n-vector a^(MOS) and an n-vector b^(MOS), the elements of a^(MOS) and b^(MOS) comprising a set of modeled output statistics (MOS) estimates of the true values of forecast correction coefficients a and b; (b) receiving, using the computer, data of an n-vector a^(prior) and an n-vector b^(prior) the elements of a^(prior) and b^(prior) comprising a set of prior estimates of forecast correction coefficients a and b; (c) receiving, using the computer, data of a spatial covariance matrix G_(a) for a^(prior) and a spatial covariance matrix G_(b) for b^(prior) (i) wherein the spatial covariance matrix G_(a) was estimated as G_(a)=E(diag(g))E^(T), where E is an n×n matrix listing a real set of discrete orthonormal basis functions for domain of interest, E^(T) is the transpose of E, and diag(g) is a diagonal matrix whose non-zero elements are given by g _(k,l) =Bexp(−b ²(k ² +l ²)) and g ^(T) =[g _(,k,l) ,k=1, . . . ,n _(k) ,l=1, . . . ,n _(l)]  in which k and l are sinusoidal basis functions with a wavenumber k in the x-direction and a wavenumber l in the y-direction, and B and b minimize the function $J = {\frac{1}{2}{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}\left( {\left( g_{k,l} \right) - \left( d_{k,l} \right)} \right)^{2}}}}$  by setting $\frac{\partial J}{\partial B} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{\frac{\partial g_{k,l}}{\partial B}\left( {\left( g_{k,l} \right) - \left( d_{k,l} \right)} \right)}}} = 0}$  and ${\frac{\partial J}{\partial\left( b^{2} \right)} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{\frac{\partial g_{k,l}}{\partial\left( b^{2} \right)}\left( {\left( g_{k,l} \right) - \left( d_{k,l} \right)} \right)}}} = 0}},$  where d_(k,l) is the element of the vector d=[E^(T)(a^(prior)−a^(MOS))]⊙[E^(T)(a^(prior)−a^(MOS))] corresponding to the wavenumber k and the wavenumber l, where the symbol ⊙ indicates the elementwise product, and where the derivatives of J are of the form $\mspace{79mu}{\frac{\partial J}{\partial B} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{{\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}\left( {{B\;{\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}} - \left( d_{k,l} \right)} \right)}}} = 0}}$      and ${\frac{\partial J}{\partial\left( b^{2} \right)} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{{- \left( {k^{2} + l^{2}} \right)}B\;{\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}\left( {{B\;{\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}} - \left( d_{k,l} \right)} \right)}}} = {\left. 0\Rightarrow{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{{- \left( {k^{2} + l^{2}} \right)}{\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}\left( {{B\;{\exp\left\lbrack {- {b^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}} - \left( d_{k,l} \right)} \right)}}} \right. = 0}}};$ and (ii) wherein the spatial covariance matrix G_(b) was estimated as G_(b)=E(diag(g^(b)))E^(T), where diag(g^(b)) is a diagonal matrix whose non-zero elements are given by g ^(b) _(k,l) =Aexp(−a ²(k ² +l ²)) and g ^(bT) =[g ^(b) _(k,l) ,k=1, . . . ,n _(k) ,l=1, . . . ,n _(l)]  where k and l are the sinusoidal basis functions with the wavenumber k and the wavenumber l, and where A and a minimize the function $J^{b} = {\frac{1}{2}{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}\left( {\left( g_{k,l}^{b} \right) - \left( d_{k,l}^{b} \right)} \right)^{2}}}}$  by setting $\frac{\partial J^{b}}{\partial A} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{\frac{\partial g_{k,l}^{b}}{\partial A}\left( {\left( g_{k,l}^{b} \right) - \left( d_{k,l}^{b} \right)} \right)}}} = 0}$  and ${\frac{\partial J^{b}}{\partial\left( a^{2} \right)} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{\frac{\partial g_{k,l}^{b}}{\partial\left( a^{2} \right)}\left( {\left( g_{k,l}^{b} \right) - \left( d_{k,l}^{b} \right)} \right)}}} = 0}},$  where d^(b) _(k,l) is an element of the vector d^(b)=[E^(T)(b^(prior)−b^(MOS))]⊙[E^(T)(b^(prior)−b^(MOS))] corresponding to the wavenumber k and the wavenumber l and the derivatives of J^(b) take the form $\mspace{79mu}{\frac{\partial J^{b}}{\partial A} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}\left( {{A^{b}{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}} - \left( d_{k,l}^{b} \right)} \right)}}} = 0}}$      and ${\frac{\partial J^{b}}{\partial\left( a^{2} \right)} = {{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{{- \left( {k^{2} + l^{2}} \right)}A\;{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}\left( {{A\;{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}} - \left( d_{k,l}^{b} \right)} \right)}}} = {\left. 0\Rightarrow{\sum\limits_{k = 1}^{n_{k}}{\sum\limits_{l = 1}^{n_{l}}{{- \left( {k^{2} + l^{2}} \right)}{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}\left( {{A\;{\exp\left\lbrack {- {a^{2}\left( {k^{2} + l^{2}} \right)}} \right\rbrack}} - \left( d_{k,l}^{b} \right)} \right)}}} \right. = 0}}};$ (d) receiving, using the computer, data of a spatial error covariance matrix R_(a) for a^(MOS) and a spatial error covariance matrix R_(b) for b^(MOS), (i) wherein R_(b) was estimated using the relation $R_{b} \approx {E\left\{ {{diag}\left\{ {E^{T}\left\{ \frac{\frac{1}{M - 1}{\sum\limits_{m = 1}^{M}{\left( {{\underset{\_}{\zeta}}^{m} - \underset{\_}{\overset{\_}{\zeta}}} \right)\left( {{\underset{\_}{\zeta}}^{m} - \underset{\_}{\overset{\_}{\zeta}}} \right)}}}{M} \right\} E} \right\}} \right\} E^{T}}$ where ξ^(m)≈f^(m)−a^(MOS)⊙y^(m)−b^(MOS) where gives a vector of differences between a forecast [(f_(m))^(T)=[f₁ ^(m), f₂ ^(m), . . . , f_(n) ^(m)]] and a MOS prediction f_(m) ^(MOS)=a^(MOS)⊙y^(m)+b^(MOS) of a forecast from an mth event and ${\underset{\_}{\overset{\_}{\zeta}} = {\frac{1}{M}{\sum\limits_{m = 1}^{M}{\underset{\_}{\zeta}}^{m}}}},$  and wherein {R_(b)}_(ij) denotes an element of the matrix R_(b) lying on its ith row and jth column; and (ii) wherein R_(a) is estimated as R_(a)≈D(Y⊙R_(b))D, where ${\left\{ Y \right\}_{ij} = \left\{ {\sum\limits_{m = 1}^{M}\frac{\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)}{M - 1}} \right\}},{{{and}\mspace{14mu} D} = \left\{ {{diag}(Y)} \right\}^{- 1}},$  and wherein each element {R_(a)}_(ij) of error covariance matrix R_(a) was estimated as $\begin{matrix} {\left\{ R_{a} \right\}_{ij} = \frac{\left\langle {\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)} \right\rangle\left\{ R_{b} \right\}_{ij}}{\left\lbrack {\frac{1}{M}{\sum\limits_{m = 1}^{M}\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)^{2}}} \right\rbrack\left\lbrack {\frac{1}{M}{\sum\limits_{m = 1}^{M}\left( {y_{j}^{M} - {\overset{\_}{y}}_{j}} \right)^{2}}} \right\rbrack}} \\ {{\approx \frac{\left\{ {\sum\limits_{m = 1}^{M}\frac{\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)\left( {y_{j}^{m} - {\overset{\_}{y}}_{j}} \right)}{M - 1}} \right\}\left\{ R_{b} \right\}_{ij}}{\left\lbrack {\frac{1}{M - 1}{\sum\limits_{m = 1}^{M}\left( {y_{i}^{m} - {\overset{\_}{y}}_{i}} \right)^{2}}} \right\rbrack\left\lbrack {\frac{1}{M - 1}{\sum\limits_{m = 1}^{M}\left( {y_{j}^{M} - {\overset{\_}{y}}_{j}} \right)^{2}}} \right\rbrack}},} \end{matrix}$  where y _(i) is a sample mean of a verifying variable y_(i) ^(m) over a plurality of m events; and (e) transforming the data of a^(MOS), b^(MOS), a^(prior), b^(prior), G_(a), G_(b), R_(a), and R_(b) into a set of estimated regression coefficients a^(FMOS) and b^(FMOS), wherein a ^(FMOS) =a ^(prior) +G _(a)(G _(a) +R _(a))⁻¹(a ^(MOS) −a ^(prior)); and b ^(FMOS) =b ^(prior) +G _(b)(G _(b) +R _(b))⁻¹(b ^(MOS) −b ^(prior)). 